The demand for safer and more efficient products is challenging engineers to reduce time-to-market, shifting prototyping from physical to virtual environments. Machine learning, particularly artificial neural networks, aids in this transition by building models that may represent complex engineering systems. In multibody system dynamics, where simulations are computationally demanding, artificial neural networks offer two key benefits: (i) surrogate modeling, providing reduction techniques for complex, non-linear systems, and (ii) handling partially unknown dynamics using experimental data. Four primary strategies have emerged for predicting time histories of motion in mechanical systems using artificial neural networks: (i) Incorporating time as an input to the artificial neural networks to predict the desired solution variables at that specific time. (ii) Utilizing an entire time series, such as loading data, as input, the network generates the corresponding time series of the desired solution variables in a single pass through the artificial neural network. (iii) Employing an artificial neural network to advance one time step into the future using the states as input. (iv) Leveraging the artificial neural networks solely for learning the equations of motion or energetic quantities, e.g., Lagrangian/Hamiltonian, coupled with standard techniques from analytical dynamics and time integration. Approaches (iii) and (iv) may be considered advantageous owing to their inherent physics-informed nature, greater flexibility in adopting varying time steps, and ease of integration with classical simulation techniques. This contribution integrates various classical time integration schemes with multibody systems modeled using neural networks, whose architecture has proven effective for mechanical systems, as shown in [Slimak et al. JCND 2024, https://doi.org/10.1115/1.4065728]. The equations of motion for these systems can be formulated in either first- or second-order forms, with neural networks approximating the right-hand side terms. The first-order form is compatible with Runge-Kutta methods, while the second-order form aligns with Newmark methods. As demonstrated in [Slimak et al. JCND 2024, https://doi.org/10.1115/1.4065728] the choice of time integration scheme is crucial for achieving accurate, long-term solutions. This work, therefore, seeks to advance the application of artificial neural networks in simulating the dynamics of complex mechanical systems.